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Fast  Evaluation  of  Radial  Basis  Functions: 

A Multivariate  Momentary  Evaluation  Scheme 

R.  K.  Beatson  and  E.  Chacko 


Abstract.  This  paper  presents  a scheme  for  fast  evaluation  of  splines, 
or  “radial”  basis  functions,  of  the  form  s(-)  = p(-)  -f  XAi  — zi)- 
Here  p is  a low  degree  polynomial  and  $ : ]Rn  — * 1R  is  a function  that 
need  not  be  radial.  This  multivariate  momentary  evaluation  scheme  is 
a generalization  of  the  fast  multipole  method  in  which  calculations  with 
far  field  expansions  are  replaced  by  calculations  involving  moments  of  the 
data.  The  primary  advantage  of  this  new  algorithm  is  that  it  is  highly 
adaptive  to  changes  in  $. 


§1.  Introduction 

This  paper  presents  a scheme  for  fast  evaluation  of  splines,  or  “radial”  basis 
functions,  of  the  form 


s(-)=P(-)+XA$('_2A  i1) 

i=l 

Here  p is  a low  degree  polynomial,  and  $ : IR"  — > It  is  a function  that  need  not 
be  radial.  This  multivariate  momentary  evaluation  scheme  is  a generalization 
of  the  fast  multipole  method  in  which  calculations  with  far  field  expansions 
are  replaced  by  calculations  involving  moments  of  the  data.  The  primary 
advantage  of  this  new  algorithm  is  that  it  is  highly  adaptive  to  changes  in  3>. 
More  precisely,  changing  to  a new  basic  function  $ only  requires  coding  a one 
or  two  line  function  for  the  (slow)  evaluation  of  $.  In  contrast,  adapting  a 
conventional  fast  multipole  code  to  a different  $ requires  considerable  analysis 
of  appropriate  expansions  and  transformation  theorems,  followed  by  writing 
a very  specific  code.  The  current  algorithm  reduces  the  incremental  cost 
of  a single  extra  evaluation  from  O(N)  to  0(1)  operations,  and  the  cost  of 
a matrix-vector  product  (that  is,  evaluation  at  all  centres)  from  £?(V2)  to 
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0(N  log  N)  operations.  The  algorithm  can  be  viewed  as  a spline  fitter  in 
that  the  approximation  it  produces  is  a piecewise  polynomial  supplemented 
by  appropriate  direct  evaluations.  The  method  described  is  a multivariate 
generalisation  of  the  method  of  [4] . 

In  outline  the  setup  phase  of  method  is  as  follows.  Firstly,  space  is  divided 
in  a hierarchical  manner.  For  example,  in  a 2-D  setting  an  initial  square 
could  be  divided  into  a quadtree.  Then  centres  are  associated  with  the  panels 
they  lie  in.  Next,  proceeding  up  the  tree  level  by  level,  the  moments  of  the 
coefficients  (the  A:’s)  about  panel  centres  are  calculated.  Next  working  down 
the  tree  for  each  level,  a number  of  approximations  to  are  formed.  Then  for 
each  panel  within  a level,  polynomial  approximations  to  that  part  of  s due  to 
far  away  centres  are  formed  by  combining  moments  and  the  approximations 
to  $.  The  evaluation  phase  first  identifies  the  childless  panel  to  which  the 
evaluation  point  belongs.  Then  it  approximates  the  far  field  by  evaluating 
the  polynomial  associated  with  that  panel,  and  adds  to  that  approximation 
the  directly  calculated  near  field  part  of  s( x).  For  reasons  of  space,  we  will 
not  detail  suitable  methods  for  subdividing  space,  or  the  process  of  evaluation. 
These  matters  are  well  understood  in  the  context  of  the  fast  multipole  method, 
see  for  example  [2], 

The  paper  is  organized  as  follows.  The  necessary  mathematics  for  form- 
ing polynomial  approximations  to  s from  moments  and  approximations  to  4>, 
and  for  translating  moments  is  given  in  Secton  3 below.  Section  4 contains 
symmetry  results  that  can  substantially  reduce  the  amount  of  work  required  to 
form  approximations  to  $ at  each  level.  Section  5 contains  numerical  results 
obtained  with  a preliminary  implementation  of  the  method. 

§2.  Notation 

We  will  need  the  following  notation.  A multi-index  a = (ai,...,a„)  is  an 
n-tuple  of  nonegative  integers.  If  x is  an  element  of  ]Rn,  we  will  write  its 
components  as  xi,X2,  ■ ■ ■ ,xn.  We  will  also  need  sequences  of  points  in  ]Rn. 
In  an  effort  to  make  the  meaning  of  all  subscripted  symbols  transparent,  we 
will  write  all  such  sequences  of  vectors  as  {zm},  and  z will  never  be  used 
unsubscripted  to  denote  a single  point  in  ]Rn. 

If  a and  6 are  elements  of  lRn,  then  we  will  say  a is  less  than  or  equal  to  6, 
and  write  a < b,  if  a;  < 6;  for  all  1 < i < n.  We  also  define  for  x 6 !Rn  and  a,  (3 

multi-indices  |a|  = aq  -ba^d h an,  a!  = aqla^!  - • • anl,  xa  = x°' x%*  • • • , 

and  for  0 < (3  < a take  = \n"pyp\  ■ The  (multivariate)  Binomial 

Theorem  then  assumes  the  form 

(*+o°  = E (oV-'v, 

0 <0<a  ' 


for  all  multi-indices  0 < a G Zn  and  points  x,t  € Hn.  We  further  define  the 
normalized  monomial 


Va(x)  = xa  / a\ 


(2) 
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M 

&t,a  “ ^ ' K,mVcc{t  Zm)‘  (3) 

m= 1 

Also  define  [a,  6]  to  be  the  n dimensional  box  {a;  £ ]Rn  : a < x <b},  and  e to 
be  the  n- vector  [1, 1, , 1]T. 

§3.  Moment  Expansions 

Lemma  1 below  shows  that  we  can  form  a polynomial  approximation  si  to 
a spline  s of  the  form  (1)  by  combining  moments  of  the  weights  A with  the 
coefficients  of  a polynomial  approximation  q to  $ on  a “double  width”  panel. 

Lemma  1.  Approximation  via  moments  - correlations.  Let  c,  d £ 1R"  with 
c,  d > 0.  Let  t £ ]Rn  with  |f,|  > c*  + d{  for  1 < i < n.  Let  e > 0 and  $ be  a 
function  in  C[t  — (c  + d),  t + (c  + d)}.  Let 

q(x)  = a<*va(x  - 1) 

{a:0<a<ke} 


be  a polynomial  of  coordinate  degree  k such  that 

11$  - 9lU°°[t-(c+d),t+(c+d)]  < e • 

Given  centres  zi,Z2,. . ■ zm  with  zm  £ \—d,  d]  C 1R"  for  1 < m < M,  and 
weights  «i,  «2,  • • • , Km  € IR,  let  the  corresponding  “radial”  basis  function 

M 

S(X)  = K™®(x  ~ zm),  (4) 

m= 1 


be  approximated  by 


M 

Sl(x)  = ^2  K™q{x  - Zm)  ■ (5) 

m— 1 


Then 

II5  S1 1 1 Z, 00 [t — c,i+c]  ^ e||^||l  • 

(6) 

Moreover, 

«i(z)  = bpVp(x-t), 

(7) 

{0-.O<0<ke} 

where 

bp  = ^ ^ aa^0,a— (3- 

(8) 

{or.  (3  < a < he} 
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Remark:  Often  in  our  applications  of  this  lemma,  will  be  truly  radial,  i.e. 
of  the  form  $(•)  = </>(||  • H2),  for  some  function  <f>  of  one  variable. 

Proof:  If  -c  < x — t < c,  then  for  1 < m < M we  have 


— ( c + d)  < (x  — zm)  — t < c + d . 


This  shows  all  but  the  expression  for  s\  in  terms  of  the  moments  of  the  data. 
To  see  the  latter,  write 


M 


«l(z)  = E 

m=l 

M 

m=l  0<a<ke 

M ( 

= E Km  { E 5 ^ " Zm)° 


m=l  ^ 0<a<ke  J 

= X>-(  £ 5 E 

m=l  I 0<a<ke  ' O<0<a 


0-/3 


JW 


= 53 "»» s 13  a“  13 

m=l  I 0<a<fce  0</?<a 


(g  - t)0  (~zm) 


CX  — fi 


f3\  (a  — /?)! 


E E Vp(x  - t)tJo,a-p  / 

0<a<fce  (o</?<a  J 

E E aa<7o,c-0  > V0(x  - t).  □ 

0 <0<ke  I 0<a<ke  J 


An  efficient  way  to  form  the  approximation  sj  is  to  use  real  FFTs  to 
compute  the  correlation  of  equation  (8).  Nominally,  such  a process  involves 
three  FFTs  per  correlation.  However,  things  may  be  arranged  so  that  the 
transforms  of  moments  of  panels,  and  those  of  the  required  approximations 
to  $,  are  computed  once  and  used  many  times.  Also  sequences  of  coefficients 
can  be  summed  in  the  transform  domain  rather  than  the  function  domain. 
This  lowers  the  average  number  of  FFTs  per  correlation  dramatically.  Thus, 
in  2 dimensions  the  work  per  correlation  is  reduced  to  approximately  0(fc2)  + 
bk2  log  k operations,  where  b is  small. 

A proof  similar  to  that  of  Lemma  1 shows  the  following  total  degree 
version  of  the  approximation  via  moments  lemma. 
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Lemma  2.  Let  c,d  € H"  with  c,  d > 0.  Let  t £ Rn  with  j<;|  > a +b,,  1 < 
i < n.  Let  e > 0 and  $ be  a function  in  C[t  — (c  + d),  t + (c  + d)].  Let 

q(x)  = ^ t), 

{a:0<a  and  |a|<k} 

be  a polynomial  of  total  degree  k such  that 

11$  - g|U~[t-(c+d),t+<c+d)]  < e • 

Given  centres  z\,  Z2,  ■ ■ . , zm  with  zrn  € [—d,  d]  for  1 < m < M and  weights 
Ki,  K2, . . . , km  € H,  let  the  corresponding  “radial”  basis  function 

M 

s{x)  = "^2  Km$(x  - zm), 

m—1 


be  approximated  by 


M 

«i(®)  = K^q{x  - zm) . 

m= 1 


Then 


Moreover, 


where 


IIs  _ Sllli^lt-c.i+c]  < e|N|l  • 

5i(s)  = J2  bpVpix-t), 

{0:O<0,  \0\ <k) 

b@  — ^ aa(7Q^a—ff  . 

{a:/3<a,|a|<fe} 


The  next  lemma  shows  that  shifted  moments  can  be  expressed  as  a convo- 
lution of  moments  about  a given  point.  This  result  will  be  used  in  generating 
the  moments  corresponding  to  a larger  panel  of  centres  indirectly  from  the 
moments  corresponding  to  subpanels.  The  indirect  process  will  be  more  effi- 
cient than  direct  formation  when  the  number  of  centres  is  large  because  the 
operation  count  for  the  indirect  shift  depends  on  the  order  of  the  moments, 
not  on  the  number  of  centres. 

Lemma  3.  Indirect  shifting  of  moments-  convolutions.  Let  z\, . . . ,zm  be  given 
points  in  ]Rn  and  ki,  . . . , km  be  corresponding  weights.  Let  oVtCl  be  the  a-th 
normalized  moment  of  the  data  defined  in  equation  (3).  Then  for  all  v,  u € 1R" 
and  multi-integers  a 


&VtCX 


= VP(V  ~ U)au,a-0  ■ 

O<0<a 


(9) 
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Proof: 


M 

“ X X V0(V  - U)V<*-0(U  - zm) 

m— 1 0</3<a 


= X VP(V  ~ U)a»,<*-P-  D 

0</3<a 


In  using  Lemma  3 to  translate  moments  we  can  reduce  the  operation 
count  by  using  either  FFT  convolution  or  a tensor  product  approach. 

We  discuss  firstly  the  tensor  product  approach.  The  tensor  product  ap- 
proach may  be  viewed  as  making  a shift  from  u to  v not  in  a single  step,  but 
as  a series  of  shifts  in  the  coordinate  directions.  For  simplicity  we  will  discuss 
only  the  2-dimensional  case. 

Consider  the  formula  (9).  If  v - u = (x,  0),  then  we  see  immediately  that 
Vp(v  — u)  is  nonzero  only  when  the  second  component  of  /?,  fo,  is  zero.  Thus 
for  v having  the  same  second  component  as  u, 


ffv,a  = X • (10) 

0</3i<«i 

Considering  aun  as  an  array  indexed  by  7 and  the  calculation  of  moments  of 
degree  not  exceeding  k,  equation  (10)  above  implies  that  each  row  of  cr„Q  may 
be  calculated  in  0(k2)  flops.  Thus  if  v and  u have  the  same  second  component 
all  moments  of  degree  not  exceeding  k can  be  shifted  in  0(k3)  flops.  Similar 
remarks  apply  to  {aw>a}  and  {crv  a}  when  w and  v differ  only  in  their  second 
components.  Thus,  in  the  2-dimensional  case  this  tensor  product  strategy 
reduces  the  flop  count  for  a single  shift  of  all  moments  of  degree  not  exceeding 
k from  0(k4)  to  0[k3). 

An  alternative  is  to  use  FFT  convolution  to  compute  the  transformation 
of  the  moment  shifting  lemma,  Lemma  3.  The  corresponding  operation  count 
is  0(k2  log  fc)  in  the  2-dimensional  case. 


§4.  Symmetry  and  Approximations  to  $ 

In  this  section  we  will  show  how  symmetry  considerations  can  greatly  reduce 
the  number  of  approximations  to  $ that  need  to  be  computed.  In  the  2- 
dimensional  situation,  with  a quad  tree  subdivision  of  space,  and  without 
clumping,  there  are  40  different  geometries  of  source  and  target  for  each  level. 
The  method  requires  approximations  to  $ on  all  the  corresponding  double 
rectangles.  However,  for  most  choices  of  $ the  number  of  approximations 
that  need  to  be  computed  from  scratch  is  reduced  to  7 by  the  symmetry 
relations  of  Lemma  4 below.  Related  symmetry  considerations  for  the  fast 
multipole  method  are  discussed  in  Wang  and  LeSar  [7]. 
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Lemma  4.  Symmetries  and  approximations  of  "radial"  functions.  Let  J\, 

Jn,  be  subsets  of  1R  and  1 < p < oo,  and  let  $ and  q be  functions  in 

• Suppose  $ is  an  even  function  of  the  k-th  component  of  x,  xk-  Define  a 


function  q by 

q(xi  j ■ • • j Xfc  > ■ • • j *^n)  — ?(*Cl  j - * • j Xfc  i • • • j *Cn ) ; (1 1 ) 

for  all  (*i, . . . ,xn)  € Ji  x . . . x — Jt  x . . . x Tien 

11$  - qWlp^  x-x-Jkx-xJn)  = 11$  - qhrtJiX-xJkX-xJn)-  (12) 

• Suppose  <b(y)  is  unchanged  by  permutation  of  the  components  j/i, . . . , yn, 
of  y.  Let  7 r be  any  permutation  of  the  integers  {1, 2, . . . , n},  let  7T”1  be 
its  inverse,  and  define  q by 

Q(yii ' ' ■ i Vn ) 4(?/ir-1(l)>  • • * > yn~1{n))i  (1^) 

for  all  y € X • • • x Tien 

11$  ~9||if(tr,(1)X-XjT(n))  = 11$  - qWLPiJxX-Xjn)'  (14) 


Proof:  Let  4>,  q and  q be  as  in  the  statement  of  the  first  part  of  the  lemma. 
Let  x = (x\,  • • • ,Xk,  ■ • • x„ ) be  a point  in  J\  x . . . X —Jk  x . . . x Jn.  Using  the 
evenness  of  $ in  the  fc-th  component  of  x, 

|#(*1,  • • • i xk  i ‘ • • ) *^n)  %n)  I 

= I$(*l,  • • • , Xfc,  • . . , Xn)  ■ • • } • • • j ^n)| 

= |$(xi, . . . , u, . . . , xn)  q^Xi, . . . , u, . . . , xn)|  , 

where  (xi,  xn)  is  a point  in  J\  x • • • X */„.  The  first  part  of  the 

lemma  follows  by  using  this  equality  in  the  appropriate  pth  power  integrals 
and  essential  supremum. 

We  turn  now  to  the  second  part  of  the  lemma.  Let  <f>,  q , and  q be  as 
in  the  statement  of  that  part  of  the  lemma.  Let  y = (yi,  ■ . ■ ,yn)  be  a point 
in  x • ■ ■ x Jn(ny  Using  that  $(p)  is  unchanged  by  permutations  of  the 
components  of  y,  and  defining  x = (xi,...,x„)  = . . . ,»»-!(„)), 

l$(2/l>  • • • > Vn)  - q{yi,  . . . , Vn) | = \&(yu  • • . ,2/n)  - q{y^(  1),  . • ■ j 2/7r-1(n) ) | 

, xn)  , • • • , #n)|  j 

where  (xi, . . . , xn)  is  a point  in  J\  x • • • x Jn.  The  result  follows  by  using  this 
equality  in  the  appropriate  pth  power  integrals  and  essential  supremum.  □ 

Remark:  Suppose  4>(-)  = <j>(\\  ■ ||p)  for  some  function  of  one  variable  <f>  and 

some  p-norm  for  H",  1 < p < oo.  Then  $ is  even  in  all  the  components  of  x, 
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Fig.  1.  Symmetry  in  the  2D  setting. 

and  is  also  unchanged  by  permutations  of  the  components  of  x.  Hence,  this 
lemma  applies  to  all  such  generalized  radial  functions  and  allows  us  to  use 
symmetry  to  obtain  approximations  to  $ on  new  regions  from  those  on  old. 

Figure  1 shows  the  geometry  of  source  panels  to  target  panels  in  1R2 
when  we  use  a quad  tree  subdivision  of  1R2  without  clumping.  The  solid 
black  square  is  the  target,  or  evaluation  panel,  and  the  possible  source  panels 
are  numbered  0 through  39.  Actually,  the  sources  in  the  left-most  column 
and  bottom-most  row  would  not  be  used  for  the  illustrated  position  of  the 
evaluation  panel  within  its  parent.  However,  they  would  be  used  for  different 
positions  of  the  target.  If  the  source  is  [— d,  d]  and  the  target  is  [t—c,  t+c],  then 
Lemma  1 requires  an  approximation  q to  $ on  the  “double  width”  rectangle 
[t  — (c  + d),  t + (c  + d)].  If  has  all  the  symmetries  of  Lemma  4,  then  at  each 
level  only  the  seven  approximations  corresponding  to  source  panels  0, . . . , 6 
need  be  calculated  directly.  The  33  other  approximations  are  easily  obtained 
by  symmetry.  The  relevant  symmetries  to  use  on  a previously  calculated 
approximation  are  indicated  in  parentheses  in  panels  7 to  39.  For  example, 
the  notation  (2*)  in  source  panel  7 indicates  that  the  approximation  for  source 
panel  7 is  obtained  from  that  for  source  panel  2 by  symmetry  in  a;.  In  the 
function  domain  this  corresponds  to  negating  the  coefficients  of  odd  powers  of 
x.  In  the  Fourier  domain  it  corresponds  to  a block  rearrangement  of  columns 
or  rows,  depending  on  which  correspond  to  x.  Similarly,  the  notation  (It) 
in  source  panel  10  indicates  that  the  approximation  for  source  panel  10  can 
be  obtained  from  that  for  source  panel  1 by  symmetry  in  x and  y.  This 
corresponds  to  a transpose  operation  on  the  coefficients  in  both  the  function 
and  Fourier  domains. 
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§5.  Numerical  Results 

Some  numerical  results  from  a primitive  implementation  of  the  algorithm  are 
given  in  Tables  1-4.  In  this  implementation  the  core  tensor  product  polyno- 
mial approximations  q to  $,  employed  as  in  Lemma  1,  are  formed  by  interpo- 
lation at  shifted  and  scaled  Chebyshev  nodes  followed  by  economisation.  The 
moments,  and  the  coefficients  of  the  approximations  to  $,  typically  have  a 
wide  dynamic  range  when  the  polynomial  degree  is  15  or  more.  Consequently, 
some  extra  device  is  needed  in  order  to  make  the  algorithm  stable,  especially 
when  the  FFT  is  used.  In  the  current  code  the  method  used  is  a scaling  of 
moments  and  polynomial  coefficients  analogous  to  that  suggested  by  Green- 
gard  and  Rokhlin  [6]  for  the  2-D  fast  multipole  method.  The  device  may  be 
viewed  as  scaling  every  panel  at  every  level  to  be  [—2, 2]2. 

In  the  calculations  reported,  the  centres  are  approximately  uniformly 
distributed  on  [0,  l]2.  If  the  number  of  centres  is  N,  then  in  Table  1 the 
number  of  evaluation  points  is  the  smallest  perfect  square  bigger  than  TV,  and 
in  Table  2 the  number  of  evaluation  points  is  the  smallest  perfect  square  bigger 
than  10TV.  In  both  cases  4>([:r,  y})  = i/c2  + x2  + y2,  where  c = 1/y/N,  and 
all  the  coefficients  Km  of  the  spline  (4)  are  1.  The  piecewise  tensor  product 
bivariate  polynomials  used  are  of  coordinate  degree  7. 


Tab.  1.  Moment  based  method  versus  direct  evaluation. 


# of 
centres 

Ord.  Alg. 
time 

FFT  Alg. 
time 

Direct 

time 

Ratio 

Abs. 

error 

Rel. 

error 

4000 

0.41 

0.27 

3.57 

13.16 

3.99E-04 

1.30E-07 

8000 

0.84 

0.58 

14.90 

25.51 

2.42E-03 

3.95E-07 

16000 

1.78 

1.23 

59.98 

48.82 

1.71E-03 

1.40E-07 

32000 

3.64 

2.65 

237.47 

89.75 

9.30E-03 

3.79E-07 

Tab.  2.  Moment  based  method  versus  direct  evaluation. 


# of 
centres 

Ord.  Alg. 
time 

FFT  Alg. 
time 

Direct 

time 

Ratio 

Abs. 

error 

Rel. 

error 

4000 

1.37 

0.97 

35.14 

36.28 

5.49E-04 

1.79E-07 

8000 

2.88 

2.00 

148.98 

74.49 

4.39E-03 

7.16E-07 

16000 

5.87 

4.17 

593.40 

142.23 

2.07E-03 

1.69E-07 

32000 

12.51 

8.98 

2356.92 

262.33 

1.82E-02 

7.43E-07 

The  timings  in  the  tables  are  in  seconds  on  an  Intel  pentium  based  ma- 
chine. Timings  are  given  for  direct  evaluation  and  for  the  algorithm  both  with 
and  without  the  speed  benefits  of  FFT  convolutions  and  correlations.  Tables  3 
and  4 repeat  the  runs  with  polynomials  of  coordinate  degree  15. 
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Tab.  3.  Moment  based  method  versus  direct  evaluation. 


# of 
centres 

Ord.  Alg. 
time 

FFT  Alg. 
time 

Direct 

time 

Ratio 

Abs. 

error 

Rel. 

error 

4000 

1.19 

0.81 

3.57 

4.42 

3.47E-08 

1.13  E- 1 1 

8000 

2.97 

1.70 

14.90 

8.74 

1.74E-07 

2.84E-11 

16000 

5.93 

3.55 

59.98 

16.88 

1.35E-07 

1.10E-11 

32000 

13.19 

7.74 

237.54 

30.69 

6.17E-07 

2.52E-11 

Tab.  4.  Moment  based  method  versus  direct  evaluation. 


# of 
centres 

Ord.  Alg. 
time 

FFT  Alg. 
time 

Direct 

time 

Ratio 

Abs. 

error 

Rel. 

error 

4000 

4.20 

2.19 

35.14 

16.06 

2.49E-08 

8.11E-12 

8000 

9.32 

4.51 

148.37 

32.89 

1.09E-07 

1.77E-11 

16000 

19.16 

9.32 

593.40 

63.65 

8.48E-08 

6.91E-12 

32000 

39.70 

19.52 

2356.92 

120.76 

3.86E-07 

1.57E-11 
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